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Abstract 

This  paper  describes  a  Monte  Carlo  sampling  plan  for  estimating  how  a  function 
varies  in  response  to  changes  in  its  arguments.  Most  notably,  the  plan  effects  this 
sensitivity  analysis  by  applying  the  acceptance-rejection  technique  to  data  sampled  at  only 
one  specified  setting  for  the  arguments,  thus  saving  considerable  computing  time  when 
compared  to  alternative  methods.  The  plan  which  applies  for  a  0-1  response  on  each 
replication  has  immediate  application  when  estimating  variation  in  system  performance 
measures  in  reliability  analysis. 

The  paper  derives  the  variances  of  the  proposed  estimators  and  allows  how  to  use 
worst  case  bounds  on  these  or  on  corresponding  coefficients  of  variation  to  choose  the 
arguments,  at  which  to  sample,  that  minimize  the  worst  case  bounds.  Individual  and 
simultaneous  confidence  intervals  are  derived  and  an  example  based  on  s-t  reliability 
illustrates  the  method.  The  paper  also  compares  the  proposed  method  and  an  alternative 
Monte  Carlo  approach  that  uses  an  importance  function. 
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Many  Monte  Carlo  sampling  experiments  aim  at  estimating  quantities  of  the  form 

g{q)  =  E  <t>(x)P{x,q)  (1) 

z€JT 

where  {0(x)}  is  a  0-1  binary  function,  (P(z,g),  is  a  probability  mass  function 

(p.m.f.)  with  given  parameter  vector  g,  and  domain  of  support  JS  so  large  as  to  make 
exact  evaluation  via  (1)  intractable.  Occasionally,  the  objective  is  to  estimate  the  function 
{g{q),  where  1=  {q  ,...,q  }.  Problems  of  this  type  arise  in  reliability  theory  where 

g{q)  represents  system  reliability  and  q.  =  (q^.^q  j  denotes  the  reliabilities  of 
components  of  types  1  through  r  which  compose  the  system  in  the  jth  of  w  component 
reliability  vectors  of  interest.  Analysis  of  giq{),---,giqw)  enables  one  to  assess  the  benefits 
of  the  alternative  reliabilities  vectors  q  ,...,qw  on  system  reliability. 

Although  one  can  simply  run  w  experiments,  sampling  from  {P(x)91)},...,{P(z,^)} 
to  produce  estimates  of  $(g  ),...,$( ^),  respectively,  a  more  efficient  method  samples  from 
{P{x,p)}  on  a  single  experiment  and  uses  these  data  together  with  the  Monte  Carlo 
importance  function  technique  or  the  acceptance- reject  ion  technique  to  produce  the  desired 
estimates.  These  approaches  are  not  new,  the  importance  function  technique  being  implicit 
in  Kahn  (1950)  and  Kahn  and  Harris  (1949)  and  the  acceptance-rejection  technique  being 
implicit  in  von  Neumann  (1949).  Beckman  and  McKay  (1987)  have  more  recently 
discussed  both  methods.  However,  until  recently  little  was  known  about  how  the  binary 
property  of  {<p(x)}  affected  the  sampling  properties  of  these  techniques  for  estimating  (1). 
Fishman  (1987)  provides  a  comprehensive  account  of  these  properties  for  the  importance 
function  approach.  The  present  paper  focuses  on  the  acceptance-rejection  method  and 
provides  a  comprehensive  description  of  the  sampling  properties  of  the  resulting  estimators 
that  exploit  the  binary  property  of  {<P(x)}  and  the  use  of  a  modified 
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p.m.f.  {C?(x,p)},  based  on  {P(z^p)}  and  information  on  bounds  for  {<p{x)}  and 
{g(q),  t0  sample  the  data.  This  last  modification  allows  the  acceptance-rejection 

method  to  work  with  considerably  improved  efficiency.  Although  the  paper  focuses  on 
applying  the  proposed  technique  to  reliability  estimation,  we  emphasize  that  the 
methodology  applies  to  the  considerably  wider  class  of  problems  with  binary  { 0(  x) } . 

Section  1  gives  basic  definitions  and  Section  2  describes  estimation  at  a  single  point. 
Section  3  then  describes  how  to  perform  function  estimation  using  the  acceptance-rejection 
method.  Section  4  shows  how  to  choose  the  design  parameter  p  to  minimize  either  the 
worst-case  variance  or  the  coefficient  of  variation  of  the  resulting  function  estimator, 
thereby  dramatically  increasing  the  efficiency  of  the  proposed  Monte  Carlo  procedure. 
Sections  5  and  6  show  that  even  in  the  worst  case,  the  proposed  technique  is  at  least  as 
good  as  crude  Monte  Carlo  sampling.  Sections  7  and  8  derive  individual  and  simultaneous 
confidence  intervals.  Section  9  illustrates  the  proposed  technique  with  an  example  and 
Section  10  compares  the  characteristics  of  the  acceptance-rejection  method  with  those  of 
the  importance  function  method. 

1.  Problem  Setting 

Consider  a  network  G  =  ( V,  $)  with  node  set  V  and  edge  set  $.  Assume  that 
nodes  function  perfectly  and  that  edges  fail  randomly  and  independently.  Let 

r  =  number  of  distinct  types  of  edges 
q  =  probability  that  an  edge  of  type  i  functions  i  =  l,...,r 

k.  =  number  of  edges  of  type  i 
*=(V"’*r) 

e..  =  jth  edge  of  type  i  j=  1  ;  i  =  l,...,r 

x  =1  if  edge  e  functions 

i)  i) 

=  0  otherwise 
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(2) 


(3) 

We  also  assume  that  G  describes  a  coherent  system.  A  system  of  components  is  coherent 
if  its  structure  function  {<(>{x)}  is  nondecreasing  in  each  argument  and  each  component  is 
relevant  (Barlow  and  Proschan  1981,  p.  6). 

Let  1  denote  a  set  of  w  component  reliability  vectors  of  interest.  Then  the 
purpose  of  analysis  is  to  estimate  the  s-t  reliability  function  { g(q ),  q  eJ}. 

2.  Estimation  at  a  Point 

Crude  Monte  Carlo  sampling  offers  a  baseline  against  which  potentially  more 
efficient  sampling  plans  can  be  compared.  Let  denote  A  independent  samples 

drawn  from  {P(x,q),  x€J>f.  Then 

ghU)  =  |  £  M.xf-'h  (4) 

is  an  unbiased  estimator  of  g(q)  with 


x.  =  E 1  x. .  =  number  of  functioning  edges  of  type  i 


‘  i=l  U 
x=  (xn,...,x  -,...;xrV...,xrk  ) 

1  r 

c5'=  set  of  all  edge  states  x 

r  k 


r  x. 


P(x,q)  =  P{xM)  =  n  n  W+(i-*J(i-«i)l  =  n  q  \i-q) 

.=i  j=  i  1  7  1=1 

=  p.m.f.  of  state  xeJf 

<f>(x )  =  1  if  the  system  functions  when  in  state  x 

=  0  otherwise 

g{q)  =  E  <fi(x)  P{xq) 
x£$ 

=  probability  that  the  system  functions. 


k~x. 

«  I 


var  ghU)  =  g(q)[l-g(q)]/K- 


(5) 


To  compute  g^q),  one  performs  K  trials  sampling  X  from  {P(z,ifc)?)}  and  evaluates 
(p(X)  on  each  trial.  The  corresponding  mean  total  computation  time  has  the  form 


T{g^q))  =  %  +  Kla^a^l  X\+az(JS;  q)] 


where 


a,(  Jf  q)  =  E  P(x,q)  C[x) 
6  lejr 


and 


C(z)  =  expected  time  to  evaluate  <t>{ x). 

The  quantities  a  a  ,  a  and  a  (J5,q)  are  machine  dependent. 

1  Z  O 

We  now  show  how  to  modify  the  sampling  plan  to  improve  statistical  efficiency 
using  information  on  bounds  as  described  in  Fishman  (1986).  Suppose  that  there  exist  0-1 
binary  functions  {<f>L{x),  xe^}  and  {<PJi x),xtJ£\  such  that 

<t>L{x)  <  <p{x)  <  (bj^x)  V  xtJS. 

Then  g(q)  has  lower  and  upper  bounds  gL{q)  and  gj^q),  respectively,  where 

9,(  Q)=  S  <t>(x)P{x,q)  it{L,U}. 

*  xa$  1 


Suppose  that  one  now  samples  independently  from  the  modified  p.m.f. 


where 


A(?)  =  gj^q)  -  gL(q). 

Then 

SK(<l)  =  3L(<l)*m^  2  «-<(,))  (7) 

1=1 

is  also  an  unbiased  estimator  of  #(g),  but  with  variance 

var  g^q)  =  {gLlq)-g{q)Mq)-gL(q)}/ 1(  <  A2(q)/4K.  (8) 

Compared  to  crude  Monte  Carlo  sampling,  one  has 

— °(9)  =  (9) 

var  gh{q)  L  ^  L  j 

>  1, 


indicating  that  g^q)  always  has  a  variance  no  larger  that  var  gh{q)- 

To  compute  g^q)  using  precomputed  bounds,  one  performs  K  trials  sampling  X 
from  {Q{x,q)}  and  evaluates  <t>{X)  on  each  trial  Here  mean  total  time  assumes  the  form 
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T{g^q))  =  *  K\fl^02\ g\+a3(JZ0Vp)/A(q)} 

where 

JtQl  =  {xeJT;  <pL{x)=0  and  ^^=1} 

and  PyPy  and  /?2  denote  machine  dependent  constants. 

Observe  that 

I\(q)  =  A'var  g^qj/vax  g^q) 

denotes  the  number  of  observations  one  would  have  to  take  with  crude  Monte  Carlo  to 
achieve  the  same  variance  that  arises  in  observations  using  {Q(x,p)}.  Then  A ^q)  = 
T(g K^{q)) / T{g ^q))  measures  the  efficiency  of  g^q)  relative  to  gh{q)  and  for  large  I< 
and  |  S\  has  the  approximate  form 

a  +a(Jf,q) /\8\  g(<l)[l-g{q)} 

A ,(*)«  - — - -  (10) 

[/?2+o3(jr0i,  q)/A(q)  \8\\  [guiq)-g(q)][giq)-gL(q)} 

a^  +  oJXq)/\S\  ] 

>  - — - -  D(q) 

[02+az(XQVq)/A(q)\g\\ 

where  (9)  defines  D{q)  >  1.  A  ratio  greater  than  unity  favors  the  alternative  sampling 
plan.  Experience  (Fishman  1986)  has  shown  this  to  be  the  case  for  moderate  and  high 
component  reliabilities  for  s-t  reliability. 
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3.  Function  Estimation  Based  on  the  Acceptance-Rejection  Method 

To  estimate  g(q)  for  each  component  reliability  vector  q  e  2,  —  {g  },  one 
can  perform  w  separate  experiments,  sampling  from  {Q(z,^)}  in  (6)  on  the  ith 
experiment  for  i  =  l,...,w.  This  procedure  incurs  the  cost  of  running  w  individual 
sampling  experiments.  However,  one  can  actually  avoid  this  cost  by  performing  a  single 
experiment,  sampling  data  from  (Q(j;p)}  and  then  using  these  same  data  to  estimate 
fiq^-^giq  )■  We  later  show  that  if  the  component  reliability  vector  p  at  which 
sampling  occurs  belongs  to  the  proposed  approach  leads  to  estimates  of  specified 
accuracy  at  a  cost  no  larger  than  that  incurred  by  performing  all  w  individual  experiments 
to  achieve  the  identical  accuracies. 

Consider  the  p.m.f. 

fix)  =  ab(x)  c(x)  xe  JS  (11) 


where 


fix)  >0,  E  c(x)  =  1,  0  <  fix)  <  1  and  a  -  1/  E  fix )  fix). 
x£$  x£Jf 

Suppose  one  samples  X  from  the  p.m.f.  {c($)}  and  Z  from  2^(0, 1).  If  Z  <  fiX),  then 
X  has  the  p.m.f.  {,/(*)}  in  (11).  This  acceptance- rejection  method  of  sampling  is  due  to 
von  Neumann  (1949).  For  the  current  problem, 


fix)  =  Q(x,p) 


b(x)  =  R(x,q,p)/ R*{q,p), 


(12) 


where 


s- 


R(x,q,p)  =  P(x,q)/P(x,p) 


n  (qjpt)  ‘ 

«=  1 


(l-q)/(l-p) 


krz. 

«  t 


(13) 


and 


R  (q,p)  =  ma xR(x,q,p)  s.t.  <pr(x)= 0  and  <pTlx)=l. 

X€<$  L  U 


The  quantity 


o  =  A(p)  fl*(M,P)/A(«) 


(14) 


denotes  the  mean  number  of  trials  required  until  one  successfully  obtains  an  X  from 

{<2(  *.?)}• 

A  small  modification  increases  the  efficiency  of  this  procedure.  Let 


R^P)  =  ma  xR(x,q,p) 
u  xe  Jr 

s.t.  <p{x)= 0  and  <p;/(x)=l 

(15a) 

R.(q,p)  =  ma  xR{x,q,p) 

1  XZ$ 

s.t.  <^(r)=0  and  <p{z)= 1 

(15b) 

F[^,q,p)  =  R{x,q,p)/Rt{q,p ) 


i  =  0,1.  (16) 


Suppose  one  samples  X  from  (Q(z,p)},  samples  Z  from  2^(0, 1)  and  determines  <fi(X). 
If  Z  <  F(X,<fi(X),q,p),  then  A  has  the  p.m.f.  {(?(*,?)}  with  mean  number  of  trials  until 


success 


9rkp)~9{p)  ,  X  9(p)-9r(p) 

“  =  A(g)  *o(w)  +  --  A(,j 


<  A(p)  R(q,p)/A{q), 


since  max  [(^p)  -  g{p),  g{p)  -  gL{p)\  <  A(p)  and  max  [i?0(?,p),  ^(ftp)]  <  R{q,p).  The 
computations  of  #0(tf,p)  and  R^{q,p)  depend  on  the  choice  of  bounding  functions 
{<t>L(x)}  and  {(pyix)}  and  axe  discussed  in  the  example  in  Section  9. 

We  next  describe  the  statistical  properties  of  data  generated  by  the 
acceptance-rejection  method. 


Theorem  1.  Let  X  and  Z  denote  samples  drawn  from  {Q(x,p)}  in  (6)  and  240,1) 
respectively.  Define  RQ  =  RQ{q,p),  R{  =  R^qj), 

Vt{x,u,q,p)  =  1  if  0  <  u  <  R(z,qyp)/R.  i  =  0,1  (17) 

»0(w,q,P)  =  9^9)  ~  A(p)flo[l-0(x)]<^(r)(z,w,?,p)  (18a) 


and 


0, (*.«.«.?)  =  S£(«)  +  vM(x,v,q,r). 


(1 8b) 


— 10— 

i.  E{[\-<t>{X)\vm{X,Z,q,p)}  =  [g^q)  -  g(q)}/A(p)R0 

ii.  =  [5(5)  -  g^/Aip)^ 

iii.  Ent(X,Z,q,p)  =  g(q)  i  =  0,1 

iv.  var  nQ{X,Z,q,p)  =  [g^q)  -  s(  «)][$(?)  -  g^q)  +  A (p)RQ\ 

=  v(<i)  +  [9^9)  -  g(q)}[&(P)R0-  A(tf)] 

V.  var  ^(X,Z,q,p)  =  [g{q)  -  ^(?)][A(p)H1  +  gL(q)  -  g{q)} 

=  K?)  +  l9i<l)  -  ^i(?)][A(p)^1-  A(fl')] 
vi-  cov[/x0(X,Z,g,p),  n^X^qj)}  =  v(q)  =  [g^q)-g{q)}[g{q)-gL(q)\. 
Proof.  Straightforwardly, 


E{[1  -^X^y^XJ^p)}  =  pr [<j>(X)  =  0,  <pQ{X,Z,q,p)  =  1] 


S 

xejr 


R{x,q,p) 

!W(I)1 


[4>^  *)-<t>L{z )] 

S(?J 


=  [9^9)  -  ^(?)]/A(p)f?0, 


P{z,p) 


establishing  i.  Part  ii  follows  in  analogous  fashion  and  the  proofs  of  parts  iii  through 
are  then  immediately  obvious. 
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Suppose  one  performs  A' independent  replications  generating  from  (6) 

and  from  ^(0,1).  Then 

K 

9iKiQ,P)  =jc  S  i  =  0,1  (19) 

>=i 

have  expectations  g(q)  with  var  gjK{q,p)  =  var  p.(X,Z,q,p)l  K.  Observe  that  the 
inequalities  var  nQ(X,Z,q,p)  >  v(q)  and  var  p^X^Z^p)]  >  v(q)  for  qt  p,  when  they 
occur,  signal  an  inflation  of  variances  over  what  obtains  if  one  were  to  sample  from 
{Q(z,q)}  directly.  Therefore,  it  is  of  interest  to  assess  how  much  these  variances  and 
corresponding  coefficients  of  variation  grow  when  using  the  proposed  acceptance-rejection 
method.  Theorems  2  and  3  provide  worst  case  upper  bounds. 

Theorem  2.  Let  X  and  Z  denote  samples  from  {Q(x,p)}  in  (6)  and  240,1) 
respectively.  Then 


[A(p)fl_l2/4  ifA(,)>A  (V)RJ2 


var  pt(X,Z,q,p)  <  M.(q,p)  = 


A(  ?)[A(p)Ai-A(g)]  if  A(?)  <  A(p)A./2 


(20) 


i  =  0,1. 

Proof.  Since  gL(q)  <  g{q)  <  g^q),  A  =  [g^q)  -  g{q)}[g{q)  -  g^q)  +  A [p)RQ]  has  its 
* 

maximum  at  g  (q)  =  gL(q)  +  max[0,A(g)  -  A(p)/ZQ/2],  from  which  (20)  follows  for  i  —  0. 
Similarly,  B  =  [g(q)  -  g L{q)][k(p) Rx  +  gL(q)  -  g(q )]  has  its  maximum  at  g*(q)  =  g^q)  - 
max[0,A(g)  -  A(Jp)/?1/2] ,  from  which  (20)  follows  for  i  =  1. 
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Observe  that  evaluation  of  (20)  for  j  =  0,1,  prior  to  sampling,  enables 
determine  which  estimator  has  the  smallest  worst  case  variance. 


Theorem  3.  Let 

7,(?,P)  =  [var  nt(X,Z,q,p)]* /[l-g(q)]  i  =  0,1. 

Then 

max  7q( Qip)  =  VQ{q,P)  =  (A(p) ]2/4[l-^  Jq)+A(p)R J[l-5r/?)] 

9L(i)<g(i)<9(j(i) 

if  ^(p)^0][l-5^?)-A(g)]  <  2A(?)[l-<^9)] 

=  A(?)[A(p)«0-A(?)]/[l-^i(?)]2  otherwise 
and 

max  7 ]{q,p)  =  V  {q,p)  =  [A(p)f? ]2/4[l-<?  (?)-A(p)fl  ][l-<7  (?)] 

^U)<?(*)<V<)  l  L  l  L 

if  A(p)^][l-^(?)+A(?)]  <  2A(?)[l-^(?)] 

=  ^(gXA^flj-At?)]/^-^?)]2  otherwise. 

o 

Proof.  We  give  the  proof  for  max  (g,p).  Let 

4  =  var  p0{X,Z,q,p)]/[l-g{q)]2. 

Then 


one  to 


(21a) 


(21b) 


(22a) 


(22b) 


(23) 


3- 


dA  _  ~g(g){2[l-<?y(g)l-*-A(p)/?Q}+2gy(g)[l-g6r(g)]-A(p)fi0[l-2g{/g)] 

[1-J(9)13 


Since  dA/dg(q)  <  0  and  dA/dg(q)  =  0  at 

9it)  =  0j/f) 


9  (?)  =  {2g(J{q)[l-gul<q))  -  A(p)fiJl-2^$)]}/{2[l-^g)]+A(p)A0]}, 


A  has  its  maximum  at  g  (q)  if  g  (q)  >  gL{q),  which  upon  substitution  of  g  {q)  for  g(q)  in 
(23)  gives  (21a).  If  g  {q)  <  gL(q)  then  the  maximum  occurs  at  gL{q),  giving  (21b).  A 
completely  analogous  result  holds  for  max  7  (q,p). 


4.  Choosing  the  Sampling  Probabilities  p 

The  results  in  Theorems  2  and  3  play  a  critical  role  in  deciding  at  which  component 

reliability  vector  p  one  should  conduct  the  Monte  Carlo  sampling  experiment.  For  each 

i=0,l,  one  procedure  finds  the  p£l  that  minimizes  max  M.(q,p)  where  (20)  defines 

1  q£&  * 

M  (q,p)  as  the  worst  case  var  p.(X. ,Z,q,p).  Then  one  uses 


P=P n 


if  max  MQ(q,pQ)  <  max  Af^pJ 
q£&  q£2> 


otherwise, 


so  that  sampling  from  {Q(x,p)}  with  p  as  in  (24)  minimizes  the  worst  case  variance  that 

2 

can  arise.  Finding  p .  takes  w  evaluations  of  M.(qyp).  Also,  note  that 


A;  =  fmin[max  MQ(q,p  ),  max  M Aq,p  )]/i>*| 
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gives  the  worst  case  sample  size  required  to  obtain  estimates  of  with 

variances  no  greater  than  a  specified  t>*.  This  valuable  information  can  assist  a  user  of  the 

Monte  Carlo  method  before  any  sampling  begins. 

The  proposed  technique  can  also  accommodate  a  relative  accuracy  specification. 

For  i  =  0,l,  an  alternative  procedure  finds  the  pE&  that  minimizes  max  N.(q,p)  where 

‘  qEl  * 

(21)  and  (22)  define  NQ(q,p)  and  N^p),  and  then  uses 

P  =  P0  if  max  N  ( q,pQ )  <  max  {q,Pl) 

q£jl  qE  H 

=  p^  otherwise. 

Sampling  from  { Q{x,p )}  with  this  p  minimizes  the  worst  case  coefficient  of  variance.  Also, 

K—  =  rmin[max  N  (q,p  ),  max  N  (q,p  )]/^]  (25) 

q$2  u  u  qE&  1 

provides  the  worst  case  sample  size  needed  to  estimate  giqj, ...,g(qw)  with  coefficients  of 
variation  no  greater  than  a  specified  u*. 

5.  Efficiency 

Naturally,  the  appeal  of  any  proposed  sampling  plan  depends  on  the  cost  saving  it 
offers,  when  achieving  a  specified  accuracy  as  compared  to  other  more  conventional 
methods.  These  cost  considerations  have  two  components,  one  based  on  variances  and  the 
other  based  on  computer  times  expended  per  replication.  Theorem  4  derives  an  expression 
for  the  smallest  variance  ratio  that  one  can  expect  to  achieve  when  comparing  a  crude 
Monte  Carlo  estimate  9hiq)  to  an  estimate  gih{q,p)  based  on  the  proposed  method.  This 
smallest  ratio  is  analogous  to  D(q)  in  (9)  and  reveals  the  least  favorable  circumstance  that 
one  can  expect  to  encounter.  The  ratio  can  be  computed  prior  to  sampling,  thereby 
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providing  a  lower  bound  on  what  to  expect. 


Theorem  4.  Let  Y  denote  a  sample  from  {P{z,p)},  X  a  sample  drawn  from  {Q(x,p)} 
and  Z  a  sample  drawn  from  #(0,1).  Let 

Bi  ( g{q),P )  =  var  <t>{  Y)/var  p.(X,Z,q,p )  i  =  0,1. 


Then 


min  BQ(g{q),p)  =  1/ {{guiq)[l-g(Jiq)+A(p)R 0]}*-{[l-<?  Jq)][gJq)-A(p)R  ]} *}5 

gL(q)<g{q)<g(j(q) 


(26) 


if  A(p)/?0  < 


^(g){^(g)[1-^{/(?)]+g{/g)[l-</L(<?)]} 


=  <7£(?)[l-i7Z;(q)]/A(?)[A(p)/20-A(?)]  otherwise 
and 


9  1  2 


min  B^(g(q),p)  =  l/{{5£(?)[l-^(?)-A(p)fi1]},-{[l-<7  (?)][<;  (<r)+A(p)^  ]}"} 


(27) 


.{A{p)R  :&(lH9L(lW-9ll(l)]*9l/.<Slt-9,(<l)}} 


=  ^[/?)[l-5i/q)]/A(?)[A(p)7?1-A(g)]  otherwise. 
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Proof.  We  prove  the  result  for  Bx{g{q),p).  Observe  that  dBJdgiq)  =  0  has  roots 


r,  =  1/j 


H(-D‘ 


1  ~9M)-^p)Ri 


gL{q)+H{p)kl 


(28) 


If  &{p)R1  <  1  - 
Since  dBJdg(q ) 


gjq),  the  roots  are  real  with  either 

Jj 

<  0,  then 

9(9)  =  9^9) 


r2  <  0  or  r2  >  1  and  gL{q)  <  r2  <  1. 


min  B^giq),?)  =  B^r^p) 

9L(9)i9(9)i9(/.9) 


if  r2<9u{q) 


=  B^g^q)^) 


if  r2>  g^q). 


Expression  (27)  follows  from  substituting  (28)  for  r2  in  the  inequality.  If  A {p)R{  >  1 
gAq),  then  the  roots  are  complex  and  dBjdg(q)  <  0  for  all  g(q)  6  [gL(q),  g^q)]  so  that  the 
minimum  occurs  at  g(q)  =  g^q).  Moreover,  complex  roots  imply  that  the  condition  in  the 
upper  branch  of  (27)  is  always  true,  thus  completing  the  proof.  An  analogous  result  holds 

for  BQ(g{q),p).  m 

The  availability  of  (26)  and  (27)  for  each  q  in  1  again  provides  valuable 
information  to  the  Monte  Carlo  user  prior  to  experimentation.  In  particular,  it  identifies 
at  which  q  adverse  variance  ratios  may  occur.  However,  measuring  the  statistical 
efficiency  of  (g^q,p),  and  {glK{q,p),  q^S}  as  estimators  of  {$(?),  qeS}  calls  for  a 
more  elaborate  analysis  than  that  for  estimation  at  a  single  point.  In  particular,  the 
sobering  observation  that  R^  and  R ^  in  (26)  and  (27)  increases  exponentially  with  |  <$\ 
makes  one  circumspect  about  the  benefit  of  the  proposed  method  as  the  size  of  G  grows. 
We  now  show  that  this  benefit  is  assured  for  finite  w  =  |  3\  and  number  of  edge  types  r, 
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provided  that  pel. 

Recall  that  1=  where  q.  =  and  is  the  reliability  assigned 

to  components  of  type  i  in  the  jth  component  reliability  vector  for  j  =  l,...,io.  Let 
and 


=  {t  €  pfq{-  for  at  least  one  j-  j=  l,...,u/}, 

so  that  j  J?*\  component  reliability  types  vary  in  1. 

Algorithm  A-R  describes  the  steps  for  computing  the  estimates  and  provides  the 

basis  for  measuring  efficiency.  In  addition  to  computing  {$0AXtf,p),  it 

computes  {V[ gQJ^q,p)\,  V[glh{q,p)\l  qtl)  as  unbiased  estimators  of  {var  5QA(?,p), 

var  glKiq,p)',  q£l\.  Observe  that  preprocessing  in  step  1  takes  0(\<&  *\w)  time, 

postprocessing  in  step  3  takes  O(w)  time  and,  on  each  replication,  sampling  in  step  2a  takes 

0(|£|)  time  using  Procedure  Q  in  Fishman  (1986),  summation  in  step  2c  takes 

0(  S  k .)  <  0(  |  8\ )  time  and  step  2d  takes  0(  |  <&*\  w )  time.  One  can  also  show  that  the 
ttcX*  ' 

mean  total  time  for  K  replications  using  Algorithm  A-R  has  the  form 

T{{gQhU,P)~9lhU,P)})  =  %  +  Wj !  <%*\ W  +  +  K 1«3  -  u4\ g\  +  az(Jf0Vp)/A{p) 

+  a;  |  J£*\  w  +  u  „  S 
5  6  «€  <%* 

time  where  u>0,...,w6  denote  machine  dependent  constants.  To  reduce  numerical  error,  all 
computation  in  step  3  should  be  performed  in  extended  precision  arithmetic. 
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Algorithm  A-R 


Purpose:  To  estimate  the  reliability  function  {$(}),  q  €.2}. 

Input:  Network  G  =  {¥,$)',  number  of  type  of  components  r  ;  it.  =  number  of  components  of 

type  i  for  i=l,...,r  ;  sampling  distribution  {  Q(x,p),  i  £  J$};  =  set  of  component  types 

that  vary  in  .2;  lower  and  upper  bounds  {g^{q)y  9 

and  number  of  independent  replications  K. 

Output:  {~gQKU,P),  J»).  ^1  M  unbiased  estimates  of 

{$(«)■  lK«).  var  90I^,q,P),  var  j^f.p);  f  €$■ 

Method: 

1.  Initialization 


a.  A(p)  e-  g^p)  -  </£(p). 

b.  For  each  qd2,. 

A(0,,)  =  K(l,q)  «-  0. 

For  each  i  €c^f*: 

«“  logt^l-P.l/p/i-flj)]  and  0.(«)  <-  log(l-g.)/(I-p.)]. 

2.  On  each  of  A  independent  trials: 

a.  Sample  X„ j  =  1,...,^  i  =  from  {(?(*, p)}. 

b.  Determine  d(X). 


c. 

d. 

e. 


* . 

For  each  i  €<#*:  X  «-  E*  X ... 


Sample  Z  from  3£(0, 1). 

For  f€J2: 

n#)  -  o. 

For  each  i  £<#*:  Tf?)  -  71(f)  +  i^(?)  +  X^f). 


A(X,f,p)  «-  exp  [71(f)]. 

FTX,d(X),,,p)  -  R{X,q,p)l RM(q,p). 
9<t>{X)(X,Z,q,p)  «-  |Z  +  F(X,d(X),?,p)J. 
*W*),f)  -  A(*(X),f)  +  ^(A)(X,Z,?,p). 

3.  Computation  of  summary  statistics 
For  each  q  €  Jl  : 


~9qRU,P)  <-  ?(/?)  “  A(p)fl0(,,p)A'(0,?)/A. 

*“  ?L(f)  +  A(p)rt1(f,p)A'(l.?)/A'. 

P[j0A4f,p)]  -  [A(p)R0(f,p)]V(0,*)/Al[I-A(0,f)/Al/(A-l). 
Vl?1AXf?)]  -  [A(p)fi1(f,p)]V(l,f)/Al[l-A'(l,f)/AV(A-I). 
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Let  us  now  compare  this  approach  to  estimating  {p(g),  q  el}  with  the  alternative 
approach  based  on  the  w  point  estimates  {<7^  q  e  1}  using  (4),  where  one  chooses  the 
sample  sizes  { H(q,p ),  q el}  to  achieve  equal  variances  under  the  two  methods.  That  is, 


var  =  «(«)[!-«( «)!/»( «J>) 


(29) 


where 


H{q,p)  =  K  \{q,p) 


and 


W__N_  ^(?)[l-5(9)] 

~  min  var  p .( X,Z,q,p )  ' 
i6{0,l}  ; 


Observe  that 


a(p>p)  =  gip)^-gip)]/[guip)-g{p)][gip)-gL{p)} 

and,  except  in  special  cases,  for  any  edge  type  i  e  <%* 


lim  \(q,p)  =  0  for  qf  p. 

k  .-*oo 

t 


Let 


A  (p)=  2  A(?,p). 
?e  1 


(30) 
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and  observe  that 


lim  A (p)  =  A (p,p).  (31) 

k  -*oo 


Therefore,  the  time  ratio 


A 


(32) 


where 


measures  the  efficiency  of  the  proposed  method  relative  to  using  crude  Monte  Carlo 

sampling  with  (4)  w  times  to  obtain  estimates  with  equal  variances  var  g  ( q )  = 

min  var  g  Aq,p)  for  each  qzS,.  As  k  increases,  (32)  assumes  the  form 

;£{0,l}  Jh 


Y.  [a2+a3(X,q)  /  ht]\{q,p) 

Aj(^,p)  s  SLSA - -  (33) 

<Va3(‘irorp)/A(^*«+u,6 


a2+a3{jr,p)/kt 

-  <VV°3<A ,.P)/A(p)*,  X(P'P) 


where  the  lower  bound  is  analogous  to  (10).  This  implies  that  one  should  expect  efficiency 
to  exceed  that  which  obtains  from  estimating  g(p)  only.  As  the  example  in  Section  9 
shows,  the  realized  efficiency  can  be  considerably  greater. 
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6.  Improving  Computational  Efficiency 

The  special,  but  common,  case  1  =  {^  <...<^}  provides  an  opportunity  for 
improving  the  computational  efficiency  of  Algorithm  A-R.  Write  F(x,i,q)  for  F(x,i,q,p ) 
defined  in  (16)  and  note  that  for  fixed  land  i 

0  <  F{x,i,qJ  <  F[x,i,qw_i)<...<F(x,i,ql)  <  1 

so  that  {F(z,i,qw_^]),  j  =  1  F(x^i,qQ)  =  1}  is  a  distribution  function  (d.f.).  Suppose 

that  one  draws  X  from  {<2(x,p)}  in  (6)  and  Z  from  ft(0,l),  and  let 

W=  min[*:  F{X,<t>{X),qw_^)  >  2].  (34) 

Then  (p^^{X,Z,q^p)  =  1  for  W<u;and  =  0  otherwise  for  j=l,...,W. 

Note  that  every  component  state  with  x.  edges  of  type  i  for  i  -  l,...,r  has  either 

{F{x,0,qw_^l)\  j=  l,...,w}  or  {fl(*l,*  j);  j=  1,. ..,«/}  as  its  d.f.  and  there  are  m  = 

2  II  *(Ar  +l)  of  these  d.fs.  If  to  is  sufficiently  small,  as  it  will  be  if  there  are  a  small 
«e  J?  ' 

number  of  component  types,  then  before  sampling  begins  one  can  compute  these  d.fs.  and 

use  them  to  create  tables  needed  for  the  cutpoint  sampling  method  (Fishman  and  Moore 

1984).  On  each  trial,  one  samples  W  from  these  specially  prepared  tables  in  0(1)  time 

* 

regardless  of  how  large  w  is.  Algorithm  A-R  shows  how  to  incorporate  this  alternative 

* 

sampling  method.  It  replaces  step  2e  of  Algorithm  A-R,  which  takes  0(/\  |  w)  time 

with  a  new  step  2e  that  takes  0(I<)  time,  thus  reducing  computing  time  per  replication. 


Algorithm  A-R 

Purpose:  To  estimate  the  reliability  function  {g{q),  q  Gj2}  where  .2  =  {^<...<9^}. 

Input:  Network  G  =  ( V,  £)  ;  number  of  type  of  components  r  ;  it  =  number  of  components  of 

type  i  for  t=l,...,r  ;  sampling  distribution  (Q(x,p),  *  €J?};  c^*  =  set  of  component  types 
that  vary  in  3>  ;  lower  and  upper  bounds  {^(9),  f  G^U{p}};  and  number  of 

independent  replications  K. 

Output:  {~gQfAq,p),  V[~gQI^q,p)],  V\~glI(lq,p)]\  qtty  “  unbiased  estimates  of 

{$(«).  iK«),  var  gQKXq,p),  var  01AX«,P);  9  6.2}. 

Method: 

1.  Initialization 

a.  A(p)  <-  g^p)  -  gL(p). 

b.  For  *=l,...,u>:  r(0,i)  =  r(l,i)  «-  0. 

2.  On  each  of  K  independent  trials: 

a.  Sample  «=l,...,r  from  {Q(x,p)}. 

b.  Determine 

c.  Sample  Z  from  2^(0, 1). 

d.  W <-  min[z  :  F{X,<t>{X),  >  2].  (Fishman  and  Moore  1984). 

e.  T(<t>(X),W)  »-  r(^(A),W0  +  1. 

3.  Computation  of  summary  statistics 

M0,f  J  «-  no,*)  and  A'(1,«J  <-  ni,*). 

For  *=0,1 : 

For  j=  K^q^^)  «-  A'0,9^)  +  r(»,w-7+l). 

For  >=l,„.,u»: 

♦*  -  &{p)RQ(q}-,p)K{0,qj)/K. 

9l( q})  +  A(p)R1(f^)A’(l>f^/«'. 

vI5ia(,/')]  “  IAW^1(*;^)1V1.«p/Al[1-^1.f/)/^/(*’-1)- 
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7.  Individual  Confidence  Intervals 
Since 


1 

where  $(•)  denotes  the  d.f.  of  the  standard  normal  distribution,  one  can  immediately 
compute  an  approximating  confidence  interval  for  g(q).  In  particular,  based  on  g  =  g^^i  q.p) 
and  Theorem  1,  one  has  the  approximating  100*(1  -  <5)  percent  confidence  interval 


1  im  pr 


l  [var  gtK(q,p)}* 


<0 


=  2*(/J)  - 


bl^^)-Mp)R^nK*P[A(p)R//^*l9^q)-^(p)R0-gJi4l*9\/K}i 

1  +  f?/K 


for  g(q)  where 


2 

0  =  (z :  —  [  e'y  /2  dy  =  1  -  6/2) 


00 


An  analogous  interval  can  be  computed  based  on  glf^q,p). 

Because  of  the  nonuniform  convergence  to  normality,  this  approach  inevitably 
incurs  an  error  of  approximation.  An  alternative  approach  avoids  this  error,  albeit  at  the 
cost  of  a  wider  interval. 

Theorem  5.  Let 


m(z,uj )  =  zlog(u>/z)  +  (1  -z)  log[(l-a/)/(l-z)j  0  <  z,  u<  1, 

let  u(z,6/2,I<)  denote  the  solution  to  m(z,u )  =  ^  log(£/2)  for  fixed  z  6(0,1]  and  £6(0,1), 


u  (z,6/2,K)  =  u(z,6/2,K) 


if  0  <  z  <  1 


otherwise. 


Then,  the  interval 


(p,uU)-&(p)R0w(l-K(0,q)/K,8/2,K),  g^q)-&(p)R0w\K(0,q)II<,6/2,K)) 


covers  g(q)  with  probability  >1  —  5 


ii.  (Ar(l,^)//r,V2,/iO,  ^i(g)+A(p)JR1o;"‘(l-/ir(l,ff)//ir,^/2,/iO) 


covers  g(q)  with  probability  >1-5. 


The  proof  exploits  the  observation  that 


P r[9fji<l)-^(p)R0  <  liQ(X,Z,q,p)  <  g^q)]  =  1 


Pr[^(tf)  ^  PxiXiZtP)  <  9LU )  +  ^(p)^]  =  1* 


The  resulting  confidence  intervals  follow  from  Theorem  1  in  Fishman  (1988). 
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Since  the  slowest  convergence  to  normality  for  g.^  q,p)  occurs  for  g(q)  close  to 
zero  and  unity  and  since  one  is  often  interested  in  g(q)  near  unity  the  wider  confidence 
intervals  that  result  from  this  approach  seem  a  reasonable  price  to  pay  to  be  free  of  the 
error  of  approximation  inherent  in  normal  intervals.  Since  (to(z,u>)}  is  concave  in  u,  one 
can  compute  the  required  roots  by  bisection. 

8.  Simultaneous  Confidence  Intervals 

Although  each  confidence  interval  in  Section  7  holds  with  probability  >  1-6,  the 

joint  confidence  intervals  for  { g(q ),  7  6.2}  hold  simultaneously  only  with  probability  > 

1  -  wS.  This  result  follows  from  a  Bonferroni  inequality.  See  Miller  (1981,  p.  8).  To 

restore  the  joint  confidence  level  to  1-6,  one  replaces  6/2  by  6/2w  in  (36a)  and  (36b)  and 

determines  the  corresponding  solutions.  The  effect  of  this  substitution  is  to  increase  the 

constant  of  proportionality  in  the  approximate  interval  widths  from  [21og(2/6)]’  to 

[21og(2w/6)]’  (see  Fishman  1988).  For  6  =  .01  and  u>=  20  one  has  [log(2u;/6)/Iog(2/6)]’  =■ 

1.25.  For  6  =  .01  and  w  =  100,  it  is  1.37  and  for  6  =  .01  and  w  —  1000  it  is  1.52.  However, 

if  1  denotes  a  continuous  region  in  the  |  S\  -dimensional  hypercube  (0,1)'^',  then  the 

resulting  confidence  intervals  have  infinite  widths  and  are  therefore  useless. 

For  the  case  Q  =  {$  <...<f  },  an  alternative  approach  derives  simultaneous 

confidence  intervals  for  {5(7),  76.2}  by  exploiting  the  fact  that  {K(0,q)/K;  j  =  l,...,u/} 

* 

and  {K(l,q^/I<;  j  =  l,...,w},  in  steps  3  of  Algorithms  A-R  and  A-R  ,  satisfy  the 
definition  of  an  empirical  distribution  function.  Since 

K  _1e/v(o,7)  =  p(o,7)  =  [g^-^q)\l^{v)RQ 


and 
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K~lEK{l,q)  =  p(l  ,9)  =  [g{q)-gL(q)  ]/A(p)«1  , 


pr{  ni[|A-(0,»i)//c-K0,«i)l  <  <y«)]}  2  >■ 


and 


pr{  n[\K(hqj)/K-p(l,q}\  <  d^S)]}  >  1-* 

where  d^6)  denotes  the  critical  value  of  the  Kolmogorov-Smirnov  distribution  for  sample 
size  K  at  significance  level  6.  Therefore, 

9^q)-^{p)RQ{qfP){K^A)IK±d^S)]  V;=l,...,w  (37a) 

cover  simultaneously  with  probability  >  1  -  6  and  similarly 

9L{l)  +  &{p)Rxi<IfP){K(]‘,<l}il  J  V  j=  1  (37b) 

cover  ^(q^,...,^?^  with  probability  >1-6.  For  6  =  .05,  lim  K^d^.05)  =  1.3581  and  for 

A'-*® 

for  6  =  .01  lim  I^d^. 01)  =  1.6276.  Since  dJ.OS)/^  .05)  <  1.013  for  I<  >  100  and 

K-ho 

dj^.01)/ d^.01)  <  1.014  for  K  >  80  (Bimbaum  1952),  little  error  arises  when  replacing 
d^.05)  by  1.3581/A^  and  d^.01)  by  1.627 6//^  above  for  K  >  100.  The  appeal  of  this 
alternative  approach  is  that  the  widths  of  the  intervals  are  all  independent  of  w.  The 
limitation  is  that  all  intervals  are  of  the  same  width.  In  practice,  one  can  compute  the 
intervals  based  on  (36a)  and  (36b)  with  6/2w  replacing  6/2  and  the  intervals  based  on 
(37a)  and  (37b),  and  choose  the  set  with  smaller  widths. 
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9.  Example 

An  analysis  of  the  network  in  Fig.  1  illustrates  the  proposed  method.  The  network 
has  30  edges  and  20  nodes.  The  example  assumes  r=l  so  that  all  edges  have  identical 
reliabilities,  allowing  us  to  write  q  =  q.  Note  that  any  other  specification  with  r  >  1  can 


Insert  Fig.  1  about  here. 


be  accommodated  easily.  The  objective  is  to  estimate  {5(9),  q  =  .80  +  .Ol(i-l) 
i  =  1,...,20}  where  g(q)  =  probability  that  nodes  s  =  1  and  t  =  20  are  connected  when  edge 
reliabilities  are  q.  For  sampling,  we  use  p  =  p,  again  merely  as  a  convenience.  The  lower 
and  upper  bounding  functions  {gL(q)}  and  {g^q)}  were  computed  beforehand  using 
edge-disjoint  minimal  s-t  cutsets  for  {gL{q)}  and  edge-disjoint  minimal  s-t  cutsets  for 
{#{/?)}>  as  in  Fishman  (1986).  To  determine  these  paths  takes  0(/|  £|)  time,  where  / 
denotes  the  size  of  the  smallest  minimal  s-t  cutset  and  to  determine  the  paths  takes 
0(|£|)  time.  The  determination  of  /ZQ  and  Rl  is  discussed  in  Fishman  (1988).  The 
evaluation  of  <j>(X)  using  a  depth-first  search  as  in  Aho,  Hopcroft  and  Ullman  (1974) 
takes  0(max(  |  £j ,  |  V\ ))  time. 

An  experiment  was  run  with  p  =  .80,  which  minimized  the  worst  case  variances  as 
in  (24),  and  with  sample  size  K  =  220  =  1048576.  Since  results  for  were 

considerably  more  favorable  than  those  for  {gXI^q,p)},  the  analysis  focuses  on  {~gQK(q,p)}. 
Table  1  shows  individual  point  estimates  and  confidence  intervals,  the  latter  having  been 
computed  as  in  (36a).  Table  2  compares  the  precomputed  worst  case  and  the  empirically 
observed  coefficients  of  variation  and  variances,  and  Table  3  shows  the  worst  case  and 
empirically  observed  variance  ratios,  where  the  variance  in  the  numerator  corresponds  to 
that  for  crude  Monte  Carlo  sampling. 

Recall  that  the  worst  case  results  can  be  computed  and  used  prior  to  sampling.  For 
example,  suppose  that  one  wants  a  coefficient  of  variation  no  larger  than  u*  =  .01  for  all 
point  estimates.  Since  the  largest  worst  case  results  in  Table  2  is  10.13,  one  would  use  (25) 
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to  compute  the  worst  case  sample  size  n**  =  1,008,016. 

Insert  Tables  1,2,  and  3  about  here. 

In  contrast  to  the  exact  results  in  col.  3  of  Table  1  which  took  slightly  more  than 
one  hour  each  to  compute,  all  results  in  cols.  4  through  8  took  74.9  minutes  to  compute  in 
total,  or  4.28  milliseconds  per  replication.  Whereas  the  calculated  exact  results  in  col.  3 
are  accurate  to  sixteen  significant  digits  (reduced  to  four  digits  here  for  comparative 
purposes),  the  confidence  intervals  suggest  an  accuracy  to  two  significant  digits  at  the  .99 
level.  If  two  significant  digits  is  acceptable  for  purposes  of  analysis,  then  the  Monte  Carlo 
approach  clearly  prevails. 

An  experiment  with  K  =  1048576  was  also  run  using  Algorithm  A-R  .  It,  of  course, 
gave  statistical  results  close  to  those  that  Algorithm  A-R  produced.  However,  it  took  36.6 
minutes  or  2.09  milliseconds  per  replication  revealing  a  substantial  increase  in  computing 
efficiency. 

10.  A  Comparison 

At  least  one  alternative  method  exists  for  using  the  data  from  a  single  experiment 
with  input  vector  p  to  generate  estimates  of  {g(q),  ?€.£}.  This  method  is  based  on  using 
the  importance  function  (13)  to  form 

=  9L(<j)  +  &(p)<t>(z)R(z,q,p) 
and 


th{x,q,p)  =  g^q)  +  A(p)[l-0(z)]R(z,g,p) 


-29- 


so  that  t  (X,q,p)  and  ipb(X,q,p)  both  have  expectation  g(q)  when  X  is  from  { Q(x,p)}. 
Fishman  (1987)  studies  these  estimates  in  detail  using  the  same  network,  and  a  comparison 
between  these  importance  function  (IF)  and  the  currently  proposed  acceptance-rejection 
(A-R)  estimators  seems  appropriate. 

For  every  qe£,  the  IF  estimators  have  smaller  variance  than  the  A-R  estimators 
do  and  both  methods  have  about  the  same  computation  time  per  replication.  If  variance  is 
the  dominant  consideration,  then  the  IF  method  prevails.  However,  there  are  other  issues 
that  also  deserve  consideration.  The  A-R  estimators  have  considerably  simpler  expressions 
for  variance  and  coefficient  of  variation  than  the  IF  estimators  do.  Also,  on  each  trial 
pQ(X,Z,q,p)  and  p  (X,Z,q,p)  for  the  A-R  approach  each  assume  binary  values  thus 
allowing  standard  techniques  of  analysis  for  binary  data  to  apply.  In  contrast  rpa(X,q,p) 

r 

and  i>.{X,q,p)  in  the  IF  approach  each  assume  0(  II  (k.+ 1)  values  precluding  the  use  of 

»=1  1 

the  simpler  analysis. 

With  regard  to  confidence  intervals,  the  A-R  approach  allows  one  to  compute 
individual  asymptotically  normal  intervals  without  nuisance  parameters  whereas  the  IF 
estimators  do  not.  For  individual  confidence  intervals  based  on  Theorem  5,  both  methods 
give  intervals  of  about  the  same  length.  This  is  a  consequence  of  ignoring  estimated 
variance  information  for  the  IF  method.  For  simultaneous  confidence  intervals  the  A-R 
method  allows  the  development  in  Section  8  when  fy—,?  >  whereas  the  IF  method  does 
not. 
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Computed  as  in  Algorithm  A-R. 
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Table  2 

Coefficients  of  Variation  and  Variances 
(. P  =  -80) 


[var  p  (X,Z,q,p)\ 

var 


q 

Worst  Case^ 

Observed^ 

Worst  Case^ 

Observed^ 

.80 

2.06 

2.02 

.207D-01 

.539D-02 

.81 

2.37 

2.35 

.262D-01 

.491D-02 

.82 

2.72 

2.72 

.322D-01 

.432D-02 

.83 

3.11 

3.11 

.383D-01 

.368D-02 

.84 

3.53 

3.53 

.435D-01 

.303D-02 

.85 

4.00 

3.93 

.455D-01 

.239D-02 

.86 

4.51 

4.46 

.443D-01 

.181D-02 

.87 

5.06 

4.95 

.406D-01 

.128D-02 

.88 

5.64 

5.44 

.349D-01 

.876D-03 

.89 

6.26 

5.93 

.282D-01 

.567D-03 

.90 

6.91 

6.37 

.212D-01 

.335D-03 

.91 

7.57 

6.74 

.147D-01 

.182D-03 

.92 

8.23 

7.08 

.928D-02 

.918D-04 

.93 

8.87 

7.37 

.520D-02 

.423D-04 

.94 

9.44 

7.41 

.250D-02 

.156D-04 

.95 

9.89 

7.29 

.982D-03 

.473D-05 

.96 

10.13 

6.76 

.287D-03 

.983D-06 

.97 

10.04 

5.95 

.527D-04 

.127D-06 

.98 

9.37 

4.59 

.415D-05 

.617D-08 

.99 

7.55 

2.69 

.396D-07 

.308D-10 

^  Computed  from  (21).  Estimated  from  data.  ^  Computed  from  (20). 
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Table  3 

Variance  Ratios 

var  g^q) 
var  90K(q,p) 


q 

Worst  Case^ 

Observed 

.80 

5.15 

6.48 

.81 

4.32 

5.88 

.82 

3.71 

5.47 

.83 

3.27 

5.20 

.84 

2.95 

5.06 

.85 

2.75 

5.07 

.86 

2.60 

5.21 

.87 

2.54 

5.61 

.88 

2.56 

6.18 

.89 

2.68 

7.06 

.90 

2.91 

8.55 

.91 

3.31 

10.95 

.92 

3.97 

14.73 

.93 

5.10 

20.81 

.94 

7.14 

34.12 

.95 

11.23 

63.00 

.96 

20.88 

149.23 

.97 

50.53 

471.61 

.98 

197.50 

2771.75 

.99 

2490.13 

67040.00 

^  Computed  from  (27).  ^  Estimated  from  data. 


